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A NEW TECHNIQUE FOR CALCULATING 
REENTRY BASE HEATING 

SUMMARY 


A theoretical analysis of the laminar base flow field of a two-dimensional reentry 
body has been formulated using Telenin’s method. The numerical method divides the 
flow domain into horizontal strips along the x-axis and represents the flow variables as 
Lagrange interpolation polynomials in the vertical coordinate. The complete Navier-Stokes 
equations are used in the near wake region, and the boundary layer equations are applied 
elsewhere. The boundary conditions consist of the flat plate thermal boundary layer in 
the forebody region and the near wake profile in the downstream region. The resulting 
two-point boundary value problem of 33 ordinary differential equations is then solved by 
the multiple shooting method using 12 segments. 

The theoretical aspects of the convergence of the present scheme are discussed 
thoroughly and are compared to the successful convergence of a smaller system; i.e. the 
two-dimensional, two-phase stagnation point flow solution. The unsatisfactory 
convergence of the present study, which is attributed to two shortcomings in the 
formulation, can be improved if the following two steps are taken. First, a variable 
transformed coordinate should be incorporated to allow different stretching in various 
segments such that the instabilities encountered can be avoided. Secondly, the Lagrange 
interpolation polynomials should be replaced by other forms of polynomials or analytic 
functions to remove the mathematical singularity at the rear stagnation point. 

The specific case considered in this report is that of vehicle reentry at zero angle 
of attack in a Mach 1 1 free stream with Reynolds number Re w ^ ranging from 

0.8 X 10 s to 1.2 X 10 5 . The base wall temperature remains constant at 255° K (46 (f R) 
and the free stream temperature is 217.43° K (392.28° R). It was assumed that heat 
conductivity and viscosity are linearly proportional to temperature, the specific heat is 
constant, and the Prandtl number is unity. The detailed flow field and thermal 
environment in the base region are presented in the form of temperature contours, Mach 
number contours, velocity vectors, pressure distributions, and heat transfer coefficients on 
the base surface. The maximum heating rate was found to be always on the centerline, 
and the two-dimensional stagnation point flow solution was adequate to estimate this 
value as long as the local Reynolds number could be obtained. 


INTRODUCTION 


With the introduction of reusable space vehicles, such as the Space Shuttle, 
minimum weight and reusability have become more important factors. To design the base 
region thermal protection system so that an undue weight penalty is not assessed to the 



vehicle, an accurate prediction of the reentry base region thermal environment is 
required. In addition for the case of the Space Shuttle orbiter, an accurate definition of 
the reentry base environment is required because the main engine nozzles are situated in 
the base region and are exposed to trapped recirculating gases during reentry. The 
purpose of this study is to provide a better understanding of the base separated flow 
region during reentry so that a more accurate reentry base thermal environment can be 
obtained. 

Atmospheric reentry involves total temperature and Mach number conditions that 
cannot be effectively simulated experimentally. Numerical schemes which can yield 
accurate solutions without requiring large storage capacity and long execution time for 
computers are desirable. One such scheme established by Telenin and Tinyakov[l] 
exploits the obvious numerical advantages of working with Cauchy-type problems for the 
present elliptic system of equations. It was first proposed for axisymmetric blunt body 
problems and later adopted for conical flow problems by Holt and Ndefo [2], It is well 
known that Cauchy’s problems are in general improperly posed for an elliptic system of 
equations. However for an a priori restricted class of solutions (such as the class of 
bounded analytic functions), Cauchy’s problems become correctly posed for the elliptic 
systems. Mathematically this means that to solve an elliptic system of equations by 
hyperbolic means would necessarily introduce the limitation that the solution can only be 
obtained in certain classes of functions, and the solution for this hyperbolic system exists 
only when the flow domain does not contain any discontinuities. Since the new 
hyperbolic system is arbitrary, in other words no characteristics exist, the integration of 
equations can be performed in any direction. Physically this is the process that allows the 
disturbances to propagate freely over the entire flow domain. 

Assume that the base region is composed of the base wall and two protruding 
shrouds (Fig. 1 ). The cavity walls, the free mixing layer, and the near wake region define 
the bounded domain wherein Telenin’s scheme applies. The Navier-Stokes and boundary 
layer equations are transformed so that the region of interest becomes a rectangle that is 
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subsequently divided into strips along the shrouds; Lagrange interpolation polynomials of 
degree four and seven are applied in the cavity and the near wake region across the strips. 
Augmented first-order ordinary differential equations are obtained. The problem is then 
reduced to a two-point boundary value problem. 

Errors committed in the arbitrary trial data increase exponentially with the 
number of trial variables and the physical dimension of the integration domain, so 
Telenin’s scheme is not immune to instability. This is especially true in the present case, 
because an almost singular layer, i.e. the base wall thermal boundary layer, exists right at 
the initial point. Because of the high flow variable gradients there, the errors introduced 
by inaccurate guessing of the initial values are amplified so rapidly that integration 
cannot be carried through this region. Such instability, which often appears in dealing 
with nonlinear problems, is the major difficulty in applying Telenin’s scheme. To handle 
this problem, the common simple shooting method is inadequate. The parallel or multiple 
shooting method proposed by Keller [3] and later developed by Bulirsch [4] is found 
effective in overcoming the instability. In essence the multiple shooting method reduces 
the integration domain length by subdividing the flow domain into a number of 
segments; each segment is treated by the simple shooting method. The guessed initial 
values are corrected iteratively by solving a linear system to satisfy the overall boundary 
conditions on both ends and to eliminate the discontinuities occurring at the segment 
junction points. 

There are several important advantages of the present scheme over a 
finite-difference type computation. It occupies one or two orders of magnitude smaller 
storage space; it consumes at least two orders of magnitude less computer time per 
iteration; the analyticity of the solution is guaranteed; and the equations are satisfied 
exactly on the strips. A simple estimate is given to support these assertions. The storage 
required for the present scheme is only that for storing variables at the intersection 
points of the strips and the segments; it is one or two orders of magnitude smaller than 
the number of grid points for the finite-difference scheme. The computation time for the 
present scheme is needed for the following three types of operations: 

1. Integration of N (number of variables) ■ S (number of strips) • M (number 
of segments) equation. 

2. Integration of M • (N ' S) 1 2 variational equations. 

3. Inversion of MN • S by N • S matrices. 

The number of operations for one iteration is then approximately M 3 N 3 S 3 

+ M(NS + N 2 S 2 * )F (number of integration step) or =10 8 with M = 12 , N = 7 , 

S = 7, and F = 10. The total computation time per iteration is about 10 2 seconds, 

while a finite-difference scheme would have to invert an MNSF by MNSF matrix or 

about M 3 N 3 S 3 F operations or 1 0 6 seconds per iteration. 
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formulation of the problem 
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where 'I'j(x) and 't'b(x) are top and bottom boundary layer edges, the region of interest 

becomes a rectangle bounded by i? = 0 , 17 = 1 , £ = 0 and the near wake region 
(Fig. 2). Replacing the first order derivatives by 



substituting these into equations (1) through (4), carrying out the transformation 
according to 
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Figure 2. Construction of strips and segments 
on transformed plane. 
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and rearranging, we obtain: 
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= (t + T??)u^ + ef 


( 10 ) 
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where subscripts \ and tj denote d/d$ and 3/drj. v(M) is the Prandtl-Meyer function. 
The subscripts b and bo indicate bottom boundary layer edge conditions at 
arbitrary $ and at £ = 0 ; similar conditions on the top boundary are denoted 
by t and to. We shall divide the domain of interest into S-l strips, as shown in 
Figure 2, and approximate the flow variables in terms of Lagrange interpolation 
polynomials across the strips; i.e. 



These expressions are substituted into equations (6) through (12) with the requirement 
that the resulting equations be satisfied identically on each line t?j . An approximating 

system of 7S first-order ordinary differential equations is then obtained for the 
approximate values uj , v i , pj , , oj , and of the dependent variables on 

the S lines; = constant. For the present case S = 7 and the rjj’s are 1, 0.95, 0.90, 
0.85, 0.73333, 0.61667, and 0.5. 
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Error Growth in the Boundary Layer 

The major difficulty encountered in carrying out the present computation was the 
instabilities. Gilinskiy and Telenin [ 5 j showed that an error caused by the approximation 
of flow variables by Lagrange interpolation polynomials across the strips may oscillate 
along the strips in the linear case. In present nonlinear systems, the error not only 
oscillates but grows rapidly to the neighboring strips. To cope with this, the author 
relied upon two fundamental tools, the boundary layer equations and the multiple 
shooting method; the latter will be discussed in the next section. 

The entire flow domain of interest was first conceived to be governed by the 
Navier-Stokes equations (6) through (14) so that the problem could be treated through a 
unified point of view. However this treatment experienced tremendous problems of 
instability because the uniform validity of the Navier-Stokes equations practically breaks 
down when dealing with a problem of extremely nonuniform grids. For high Reynolds 
number flows with a large separation bubble, the boundary layer equations are more 
feasible for the high gradient areas. Although this will limit the accuracy of the solution 
to less than 1/Re there, the instability problem can be avoided partially. 

The problem of error propagation in the base wall boundary layer is of special 
importance because almost all the physical processes in determining the base flow heat 
transfer properties occur there and in the free shear layers. The governing equations for 
flow can also be regarded as the error propagation equations, since without knowing the 
solution a priori the guessed initial values may contain an error of their own magnitude. 
We shall focus our attention upon the error growth of the heating rate across the base 
wall boundary layer. The following equation gives the growth rate of j3 OC 3T/dx at 
1 = 0 along the strip: 


= y (T + - 0 2 - (y - I)M®, Pr T a 2 


+ Re Pr F (p ,T , e , a ,0) (16) 

The last term on the right side can be neglected if boundary layer equations are used; 
however, if it is retained on the Navier-Stokes equations, rapid error amplification caused 
by this term will occur since the initial values cannot always be chosen so as to guarantee 
the last term’s smallness. The second and third terms are dominant then and remain to be 
negative, this will therefore reduce the danger of divergence. The approximate solution of 
the above equation can be represented by the following relation: 


2 _ 1 

p V (7 - 1 )M^, Pr T a 2 ^ V (7 - 1 )M^ Pr T a 2 ~ 
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so that (3 will decrease when £ increases; 
in other words the integration is stable. 
Similar analyses can also establish the fact 
that in shear layers the error is amplified 
slower by boundary layer equations than 
by Navier-Stokes equations. Based upon 
this result, we shall use boundary layer 
equations on base wall, shrouds, and in 
free shear layers and Navier-Stokes 
equations in the remaining regions. This 
mathematical model is depicted in 
Figure 3. 



Figure 3. Governing equations and segmen 
tation of the base flow regions. 


Boundary Layer Equations 

In the base wall thermal boundary layer, equations (11) and (12) are 
supplemented by the following governing equations: 
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( 22 ) 


pufu^ = (r + nfHpu^ + p^u) - ufp| - (p^v + pv^) 

In the forebody thermal boundary layer and shear layer downstream of the shrouds 
the following set of equations is applied: 

puf 2 U£ = puf(r + r??)^ - pvu^f + pjUjf 2 U]^ 

+ Re ( Vu + Tl W > ' (23) 
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O 

II 
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Substituting relation (15) into equations (18) through (22) and equations (23) 
through (26), we obtain an approximating system of 7 X 4 first order ordinary 
differential equations for the base wall thermal boundary layer, 7X4 equations for the 
forebody boundary layer, and 7S equations for the shear layer. Along the 
strip 7 ? - rj^ in segment 3, where the Navier-Stokes equations are applied on and below 
it and the boundary layer equations are applied above it, the vertical derivatives are 
calculated by using the same seventh order Lagrange interpolation polynomials so that 
the vertical derivatives are continuous across this strip. In the first two segments, the 
vertical derivatives are computed separately by two fourth order Lagrange interpolation 
polynomials in the forebody boundary layer and in the flow underneath the shroud. 

With the present formulation, all the following physical phenomena have been 
taken into consideration: (1) the interaction between the inviscid flow and the viscous 
flow is defined by the free interaction equations (13) and (14) along the external edge of 
the viscous layers; (2) the interaction between the shear layer and the recirculating core is 
accounted for by enforcing the continuity of the flow variables and their vertical 
derivatives 3/3ij across the strip in segment 3; (3) the upstream propagation of pressure 
wave through the shear layer near the trailing edge of the shroud is implicitly included 
through the application of Navier-Stokes equations in the near wake region (segments 2 
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and 3) and the iterative numerical scheme ; and (4) the existence of the base wall thermal 
boundary layer is explicitly formulated using the boundary layer equations along the base 
wall. If the validity of these equations is questioned near the upper left corner 
underneath the shroud in segment 1 , because of the nature ol the boundary conditions 
applied there, either the Navier-Stokes or the boundary layer equations would yield 
approximately the same results. 

To employ the present scheme, we should remember that no singularity can be 
allowed in the domain of interest except right at the segment junction points where the 
Poincarb analysis can be carried out in advance. From equations (23) and (24) it is clear 
that the rear stagnation point is a singular point of the differential equations in segments 
4 though 11. Since the location of the rear stagnation point is not known a priori , it 
poses a serious problem, because during the iteration this point may emerge in the 
integration domain such that the integral curves near it contain errors of an unacceptable 
magnitude. This difficulty can be avoided if the Lagrange interpolation polynomials are 
replaced by other forms of analytic functions; however the advantage of having the 
ordinary differential equations written in explicit form is lost. We will pursue only the 
solutions upstream of the rear stagnation point. 


Boundary Conditions 

In itial Thermal Boundary Layers . External to the shrouds, boundary conditions 
correspond to the solutions of forebody thermal boundary layers. Assuming no separation 
is ahead of the shroud trailing edge, the solutions of the compressible boundary layer 
past a flat plate are applied. With Prandtl number equal to unity, we have 



Near Wake. Solutions. The downstream boundary condition is a near wake profile 
obtained by extrapolating Kubota’s* far wake solution upstream. Defining x , y as the 
coordinates after a Stewartson-lllingworth transformation and with the origin set at the 
neck, we have 


u 

u, a 


A 



* Kubota, T.: Laminar Wake With Streamwise Pressure Gradient, GALCIT Hypersonic 
Research Project. Internal Memorandum No. 9, May 1962. 
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and 0 is the momentum thickness. Since the external flow is represented by the 
Prandtl-Meyer solution in the present study, the neck condition corresponds to that when 
the flow is parallel to the centerline. For the present problem, with M m = 11 , 

thC Me at neck = 9 586 ‘ The t0tal heat 1oSS of the flow past the vehicJe is estimated by 

neglecting the base heat transfer and assuming the vehicle length is 20 times the base 
height, so that the Stanton number is taken to be -0.0856155. The (0/H) at neck was 

put equal to 0.0618602 and ^/H = 0.6389918. Figure 4 shows the near wake profiles. 


NUMERICAL PROCEDURES - MULTIPLE SHOOTING METHOD 
AND CONTINUATION METHOD 


A serious shortcoming of the shooting method becomes apparent when the 
differential equations amplify the errors so rapidly that divergence occurs before the 
initial value problem can be completely integrated. This may happen even though 
accurate guesses are made for the initial values. The multiple shooting method can 
frequently circumvent the difficulty, or else a finite difference scheme can be employed. 
The method is essentially a combination ol difference scheme and initial value problems. 
It is designed to suppress the growth of the errors in the trial integral curves by dividing 
the domain of integration into a number of subintervals, integrating each individual initial 
value problem over its own interval, and then simultaneously adjusting all the guessed 
initial data to satisfy the boundary conditions and continuity conditions at the junction 
points. 


The formulation of the multiple shooting method can be found in Osborne [6] , 
and a comprehensive version was given by Bulirsch[4], For completeness and 

convenience in discussing the continuation method later, it will be mentioned briefly 
here. 
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Figure 4. Near wake boundary profiles. 


For a given boundary value problem, 


f - f« , y) <»> 


and 

r[y(a),y(b)l = 0 (30) 

with a < % < b , and y and r are vectors of NS components. We divide the domain of 
interest a < £ < b into M - 1 subintervals, guess a set of initial values Yj = 

every interval j = 1 - 1, and we have a set of initial value problems 
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(31) 


d$ 


f(S ,y) 


y(fj) = Yj , j = l, ... ,M-1 . (32) 

In general the function values y(?j+j , Yj) resulting from solving this system of equations 
do not match with the assumed Yj + j values at £j + j , so there exists the jump 

hj = yft j+1 , Yj) - Y j+1 , j = 1 , 2 , ... , M - 2 (33) 


and 


h M-l " r t Y i .y(?M > Y M-P ] - ( 34 > 

If we can find a set of initial values such that all hj’s and h^j. j’s are equal to zero, then 
the problem is solved. To find this solution consider the set of equations hj as functions 
of Yj and Yj + j , hj(Yj , Yj + j) # 0 . We seek a value for AYj such that hj(Yj + AYj , 
Yj + j + AYj + j) = 0 . Taking a Taylor’s expansion around Yj and Yj + j , we have 

3h { 

hj(Yj,Y j+1 ) + GjAYj + 3Y^-AY j+ , - 0 , j=l,...,M-2 


and 


'M 


-1< Y 1 >y M -i> + §1 ay, + 


3r 


3Y 


AY 


M-l 


M-l 


where 


(35) 


_ ay«j+i , Yj) 

°j sy; 
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aY j+, 


- 1 


and 


dr 

9Y M-1 


dr 9 y (t M ’ Y M-1 } 
9 Vb aY M-l 


Defining A - , B = ^ 


the correction vector A Y 


J 


we have the following linear system to solve for 



-I 

G2 "I 

0 


g M-2 



This matrix consists of M - 1 blocks and each block is of the order NS ; the order of the 
system is NS X (M - 1) . It is noteworthy that Bulirsch {4] showed that this system can 
be reduced to a set of NS X NS linear systems. Multiplying the jth block 
by . . . Gj + | for j = 1 , . . . , M - 2 and the M - 1 block by I and adding, we 

have for the unknown vector AY, , 


EAY, = R , (37) 

where E = A + BG m .j . . . Gj and -R = h M .j + BGj^.j h M _ 2 + . . . + BG m _j . . . 
G 2 G, . This linear system is solved by Gauss-Jordan elimination and the AYj’s are 
obtained subsequently by simple matrix multiplications; 

AY j+1 = hj + GjAYj , j=l, ... ,M-I . (38) 
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Common modified Newton’s method computes the initial values Yj(® + ^ 
after £ iterations by 

Yj ( e+ D = Yj® + X^AYjW , j=l , ... , M-l , and 0 « X (C) < 1 , 

however, constant X^ for all Yjj’s , i = 1 , . . . , NS would only allow a few variables 
to change while refraining the rest from varying. To accelerate the convergence, we found 
that a diagonal matrix Xjj^ was more effective here. Bulirsch (4] gave a detailed 
description of numerical computations of the matrix Gj’s and the application of 

Broyden’s technique [7] ; these will not be iterated here. The following brief discussion of 
the convergence of the shooting method given by Meng [8] however will be included for 
completeness. 

Let r and y be the boundary conditions and the unknown vector; therefore, the 
convergence sphere and rate of convergence for the shooting method are 


(1 - V' 1 - 2h 0 )7 ?0 h 0 


and 


(2h 0 ) 2c “ 1 no/2 8-1 
with the Jacobian matrix G, 


II G II * 

6 B o 

NS 

y 

9 2r i 

LS=1 

9yj 9y s 


for all i’s, 

. II G' 1 r II < rj 0 
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II r II = max I r I 
l«i«NS 


NS 

II G II = max ^ 1 G ik 1 > 

Ki<NS k~l 


and 8 is the number of iterations counted after the initial values fall within the 
convergence sphere. By the Kantorovich theorem [9] , the convergence is guaranteed as 
long as h 0 = B 0 r? Q K is smaller or equal to one-half. For simple problems, convergence 

can often be obtained by simply going through many iterations. In complex problems, 
one has to modify the guessed values to fulfill as many of the Kantorovich sufficient 
conditions as possible for convergence. Ironically, the labor required.to make such a test 
is NS times more than that needed for solving the problem itself. For example, the 
quantity K needs integration of 
MNS(N 2 S 2 + 2NS-1 )/2 equations 0 6 

throughout the entire domain so that the ^ 
advantage of working with the Cauchy-type ° 
problem will be greatly diminished. Since g 
the Kantorovich h c cannot be obtained -• 

economically for the present problem, to ^ ^ 

illustrate how the multiple shooting method < 
converges according to the theorem, we “ -3 

carried out a two-phase stagnation point oc 
flow solution. This was a smaller system of uj 
seven equations and tour subintervals; the 0 

Euclidean error norm and h„ are presented o’ 

o 

in Figure 5. One finds that the method does - 1 -9 
converge. Even the First guess falls outside 
the convergence sphere; as soon as it hits 
inside the sphere, the convergence is 
reached. Figure 5. Error norm and Kantorovich h Q . 

The subdivision of the domain for the multiple shooting method is determined by 
the relation |£j +1 - |j| ~ jj , if all the derivatives are Lipschitz continuous so that a 

stable integration can be guaranteed. In selecting this Lipschitz constant Lj , it is clear 
that the maximum row in the matrix Gj , 



1 2 3 4 5 6 7 

NUMBER OF ITERATIONS 



max 

l«i<N 
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is a measure of maximum local error growth resulting from small perturbations at the 
initial point of the subinterval j . Direct substitution of this value of B Q . to determine 

the domain length, i.e. A£j ~g , however does not yield a practical answer for the 

... °j 

multiple shooting method. Because the value of B q is quite large in nonlinear problems, 

tor example, it is of the order of 10 6 for the base wall boundary layer and of 10 3 for 
the downstream regions. Therefore in theory, about 10 3 or 10 6 subintervals to insure 
against the instability are required, but in practice the advantage of the Cauchy-type 
problem will be offset if the number of subintervals becomes comparable to the number 
of the grids by the difference scheme. This dilemma can be resolved by incorporating the 
continuation method developed by Roberts and Shipman [10] with the multiple shooting 
method. They employed the simple shooting method and stretched the domain length to 
the final length in each iteration to solve a problem which could not be solved by the 
shooting method alone. It was shown [9] that the method will be stable if the stretched 

length is bounded by > M is the uniform bound of the derivatives 

over [ Sj’Sj+l ] new and 


K = max V 

a Gj 

N 

= max ) 


l«i«N u 

9Y js 

Lt 

Ki<N k,s=l 

9Y js 9Y jk 


However it is found that one should not continue the segment length this way in practice 
either because the denominator is very large, ~0(10 10 ), but should find the Afj new by 

A ^jnew = A ^jold^ MKB Oj- > old^ MKB Oj^new > once one can have a stable integration over 

Ihe 4{ 0 | d . By the present experience, 4| jnew - « w was found 

adequate in stretching the domain length during each iteration. 


In summary by applying Broyden’s correction technique, the convergence factor, 
and the continuation method to the multiple shooting method, the present problem was 
solved using 12 subintervals during the first few iterations. In following iterations closer 
to the convergence, the number of segments was reduced to eight without any effect 
upon the stability. 

Finally it should be noted that the success of the present iteration scheme relies 
leavily upon the accuracy of the integration routine; a seventh order Runge-Kutta 
scheme with stepsize control established by Fehlberg[ll] was hence applied in the 
present analysis. 
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RESULTS AND DISCUSSION 


The present problem is reduced to a system of 33 equations after applying the 
symmetry condition along the centerline and the interaction equation along the shear 
layer. The computation was conducted on a UNIVAC 1108 computer. The program 
occupies a 54K storage space. In initial trials six segments were employed, and the 
convergence appeared poor. Later double precision and 12 segments were used, this 
improved the convergence. The bulk ot the computation time was spent in generating the 
Jacobian matrix, nearly 5 minutes each time; however by employing Broyden s technique, 
the time was reduced to 12 minutes to complete four iterations. The Jacobian matrix was 
first computed every five iterations with Broyden’s technique applied accordingly; the 
solution yielded obvious errors. It appeared that the method produces the best result if 
the Jacobian matrix is computed every three iterations. 

M 

The Euclidean error norm ^ I hj | 2 and the variation of the smallest element in 
the diagonal matrix are shown in Figure 6. The error decreases steadily for the first 
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] 3 iterations then oscillates, and the convergence factor shows similar features. This lack 
of convergence was conjectured primarily because of fixing the segment lengths by the 
considerations of stability alone, as was outlined previously. Since each set of differential 
equations applied in different segments has its own physical capabilities or limitations to 
generate certain flow patterns, the associated segment lengths over which these equations 
are integrated would therefore play an important role in achieving convergence. Based 
upon this the first three segment lengths were then treated as additional unknowns; hence 
the system was augmented to 36 equations, and the problem was solved as a multiple 
free broundary-value problem. The convergence is improved (Fig. 6), although the 
oscillation still persists. The final segment configuration indicates that the base wall 
boundary layer thickness equals 0.07478H, the two Navier-Stokes equation segments are 
of 0.004784H and 0.01 11 8H respectively followed by eight boundary layer equation 
segments of 0.1 6H each. The segment lengths vary little with the free stream Reynolds 
number. The continuation method which was mentioned in the last section succeeded in 
stretching the whole domain length from 0.2H to 0.68H smoothly. 

The evolution ot the velocity vector and the temperature contour through 
iterations to satisfy the boundary conditions and continuity across the subintervals are 
shown in Figures 7 through 9 for a Reynolds number of 10 5 . The initial velocity vector 
plot ol the forebody boundary layer revealed no turning action around the edge of the 
shroud, because the pressure drop had not propagated upstream. It is seen that there are 
discontinuities across the intervals and significant interpolation errors along 
y/H - 0.57 and 0.71 . After several iterations, this error diminished in magnitude while 
the flow around the shroud edge began to incline towards the wall. The recirculating flow 
pattern is shown clearly in Figure 9. In the initial temperature contour plot, there are 
negative values of temperatures which are indicated by the blanks, and discontinuities 
also exist across the intervals. The fact that the contour lines failed to be normal to the 
centerline indicates that errors resulting from the Lagrange interpolation exist. In an area 
near the wall, the gap between the contours is small because of the high heating rates. 
The cold and hot spots emerge in Figure 9 and the profiles show little variations along 
the horizontal direction throughout the near wake region. The diamond shape of the 
temperature contour in Figure 9 shows that some discontinuities exist across the segment 
even though the rest of the region shows quite good results. From the temperature 
contour plots, it is clearly seen that because of the small flow velocities in the base wall 
boundary layer, the heat is transferred almost entirely by the conduction process and the 
flow convects the heat generated in the shear layer to the compression region and 
recirculates it back along the centerline. The effects of the protruding shrouds upon the 
base thermal environment would be to pull the pressure rise and high heat flux occurring 
in the recompression region away from the base wall so that the heating problems to the 
base wall and engines are reduced. 

Since the plot routines pick up values only at equal vertical intervals, the flow 
variables on the shear layer edge are often missed in the temperature contours and 
velocity vector plots. The edge Mach number and pressure distributions are given in 
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Figure 9a. Velocity vector in the two-dimensional Space Shuttle base region (final iteration). 




CONTOUR IDE 

1 

.250 

2 

.750 

3 

1.250 

4 

1.750 

5 

2.250 

6 

2.750 

7 

3.250 

8 

3.750 

9 

4.250 

A 

4.750 

B 

5.250 

C 

5.750 

D 

6.250 

E 

6.750 

F 

7.250 

Q 

7.750 

H 

8.250 

1 

8.750 

J 

9.250 

K 

9.750 

L 

10.250 

M 

10.750 

N 

11.250 

0 

11.750 

P 

12.250 

Q 

12.750 

R 

13.250 

S 

13.750 

T 

14.250 

U 

14.750 

V 

15.250 

w 

15.750 

X 

16.250 

Y 

16.750 

z 

17.250 




Figure 10 for Reynolds number Re^ - 10 s - The pressure first drops smoothly along 

the forebody boundary layer and then reaches the base pressure drastically near the 
trailing edge of the shroud, while the edge Mach number increases in the same manner. 
The pressure distribution on the centerline is also shown in Figure 10; it is nearly 
constant throughout the cavity region until near x/H = 0.1 where it begins to follow 
the external pressure very closely. The Mach number at the neck, =9.587 , is also 
marked in Figure 10; if the external flow were parallel to the centerline, the external 
Mach number should be equal to this value. The vertical pressure distributions at various 
axial locations are shown in Figure 11. The pressure in the forebody boundary layer and 
downstream shear layers is nearly constant except for interpolation errors, the pressure in 
the cavity on the base wall is nearly four times higher than that at the trailing edge. 
At x/H = 0.1, the pressure drops drastically underneath the shroud trailing edge and the 
value on the centerline is close to that on the external edge. 

The heat transfer coefficients based upon the recovery temperature are shown in 
Figure 12. The maximum heating rate is always on the centerline and its value increases 
monotonically with the free stream Reynolds number. Detailed flow patterns for various 
Reynolds numbers are given in Figures 13, 14, and 15. From the velocity vector plots of 
Figures 13a, 14a, and 15a, we can find that the shroud edge Mach number varies slightly 
and monotonically with the free stream Reynolds numbers and the vector plots are quite 
similar even though the convergence was poorer for higher Reynolds numbers. The 
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Figure 1 0. Mach number and pressure distribution 
along the shear layer edge. 
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Figure 1 1 . Vertical pressure profiles at various axial locations. 



Figure 1 2. Heat transfer coefficient 
on base wall versus Y/H. 


temperature contours show that the dia- 
mond shape discontinuity disappears for 
the lower Reynolds number case. The 
Mach number contour is also shown for 
the Re^ ^ = 0.87 X 10 s case in Figure 

13c. It is seen that the major portion of 
the recirculating core is subsonic, the 
sonic line extends from the wake into the 
forebody boundary layer, and external to 
it the viscous layer is entirely supersonic. 
The subsonic region in the forebody 
boundary layer is very thin so that only 
few upstream propagating waves can 
be transmitted through this viscous 
layer. This explains the weak upstream 
influence observed in the experimental 
study of near wakes by Batt and Kubota 
[12]. The fact that a significant portion 
of the viscous layer is supersonic also 
confirms that the imbedded shocks will 
emerge deeply in the viscous region as was 
suggested by Weiss and Weinbaum [13]. 
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Figure 14b. Temperature contour in the two-dimensional Space Shuttle base region for Re - 1.1 X 10 s . 
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Figure 15a. Velocity vector in the two-dimensional Space Shuttle base region for Re - 1.2 X 10 s . 







Although we do not intend to dwell on the aspects of comparing the present 
results to experimental data since the latter do not exist for the present flow conditions 
and geometry, it is interesting to note that the present results are an order of magnitude 
different from those of Larson et al. [14]. They found Nu 2 100 ~ 150 for 
T w ~ 0-34 ~ 0.716 at = 3 and Re^^ ~ 10 7 , while the present study gives 

Nu = 5 for T w = 0.046 T^ at M«, = 11 and Re M = 10 s . The observed trend that 

the base wall thermal boundary layer thickness varies with both the Reynolds number 
and temperature difference between the wall and recirculation region is believed to be 
correct. Figure 1 6 shows the value 


Nu 

x/Re 


( T aw - T w ) 


0T, 



Re 


«,H 


obtained in comparison with the similar solution for a two-dimensional stagnation point 
flow solution given by Cohen and Reshotko [15] . The values scattered around the 
theoretical value «0.506 , and they reveal no strong dependence upon the local edge 
Reynolds number. The two-dimensional stagnation point flow is therefore seen as a close 
approximation of the base flow. Furthermore as shown by v/U^ versus y/H on the base 
wall boundary layer edge in Figure 17, the magnitude of the vertical velocity decreases as 
the Reynolds number increases, and the linear dependence upon the coordinate y/H is 
true only near the centerline. For practical purposes, it can be asserted that the base flow 
is a stagnation point type flow; the maximum heating rate can be derived from the 
stagnation point flow results so long as the local Reynolds number can be estimated. 
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Figure 1 6. Comparison with the two-dimensional 
stagnation point flow solution. 




For future investigations of even 
higher Reynolds number flows by the 
present method, the following consider- 
ations should be noted. First the insta- 
bilities encountered for integration across 
the base wall thermal boundary layer can 
be avoided if a variable transformed coor- 
dinate is incorporated to allow different 
stretching in various segments. Secondly it 
should again be emphasized that when 
using Lagrange interpolation polynomials 
in formulating a Cauchy problem, there 
should exist no singularity in the flow 
domain, because when such singularity 
.05 0 .05 .10 emerges, the advantage of using the 

v/ujd Lagrange interpolation polynomials will 

be lost. Replacing the Lagrange inter- 
Figure 17. Velocity (v) near the edge. polation polynomials by other sets of 

of the base wall boundary layer. polynomials or analytic functions can 

remove this singularity, but one more 
matrix inversion to obtain the system of 
first order differential equations would then be necessary. To include the solution 
downstream of the rear stagnation point, such replacement is needed. 

Although for simplicity we have concentrated on the zero angle of attack case, 
extensions to skew cases offer no difficulty. The various aspects of the rate ot 
convergence, the storage requirement, the computation time, and the exactness of the 
solution on the strips should cause one to favor the present method over many existing 
schemes in dealing with high Reynolds number flows. 

George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama, December 30, 1972 
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